//This function calculates the energy associated with a change in volume
void NPT() { 

	prec rho, avgrho, randno;
	prec acceptRatio, xdl, deltaL, pot, pot_old, enthalpy;
	int moveSelect=0;
	int ref = 0;
	pot_old = potentialEnergy; 

	moveCountNPT++; 

	//calculate new boxlength
	randno = genrand_real();
	deltaL = deltaLmax*(0.5-randno);
	  
	box_old.x = box.x;//save old boxlength
	box_old.y = box.y;//save old boxlength
	box_old.z = box.z;//save old boxlength

	box.x = box.x + deltaL;
	box.y = box.y + deltaL;
	box.z = box.z + deltaL;

	//calculate volumes with new and old Bl
	prec Vold = (box_old.x*box_old.y*box_old.z)*VOLUME_CF;//vol in m^3a
	prec Vnew = (box.x*box.y*box.z)*VOLUME_CF;//vol in m^3
	
	//save old positions and calculate dCM
	for (int i = 0; i<=(nMolecules-1); i++) {
		for(int j = 0; j<=2; j++) {
			//save old positions
			refConfig[i*3+j].x = hAtom[i*3+j].x;
			refConfig[i*3+j].y = hAtom[i*3+j].y;
			refConfig[i*3+j].z = hAtom[i*3+j].z;
		}
	}

	for (int i = 0; i<=(nMolecules-1); i++) {
		for(int j = 0; j<=2; j++) {
			//get distance ref. to COM
			dCM[i*3+j].x = hAtom[i*3+j].x - hAtom[i*3+2].x;
			dCM[i*3+j].y = hAtom[i*3+j].y - hAtom[i*3+2].y;
			dCM[i*3+j].z = hAtom[i*3+j].z - hAtom[i*3+2].z;
		}
	}
	//scale COM to new boxlength
	for (int i = 0; i<=(nMolecules-1); i++) {
		hAtom[i*3+2].x *= (box.x/box_old.x);
		hAtom[i*3+2].y *= (box.y/box_old.y);
		hAtom[i*3+2].z *= (box.z/box_old.z);
	}  
	for (int i = 0; i<=(nMolecules-1); i++) {
		for(int j = 0; j<=2; j++) {
			//get scaled positions by adding dCM to scaled COM
			hAtom[i*3+j].x = dCM[i*3+j].x + hAtom[i*3+2].x;
			hAtom[i*3+j].y = dCM[i*3+j].y + hAtom[i*3+2].y;
			hAtom[i*3+j].z = dCM[i*3+j].z + hAtom[i*3+2].z;
		}
	}

	//Place everything back into the box in case it moved out
	normalizePositions();
      
	//calculate new energy
	pot = calculatePotential(ref, moveSelect, hAtom); 
	enthalpy = atomPressure*(Vnew-Vold)+(pot -pot_old)-(double)nMolecules*tempK*BOLTZMANN_C*log(Vnew/Vold);

	// Accept or reject 
	randno = genrand_real();
	if(exp(-1*enthalpy/(tempK*BOLTZMANN_C))>randno) {
		acceptCountNPT = acceptCountNPT +1;
		potentialEnergy = pot;
	} else {
		Vnew = Vold;
		box.x = box_old.x;
		box.y = box_old.y;
		box.z = box_old.z;
			    
		//scale configuration ot old box size
		for (int i = 0; i<=(nMolecules-1); i++) {
			for (int j = 0; j<=2; j++) {
				hAtom[i*3+j].x = refConfig[i*3+j].x;
				hAtom[i*3+j].y = refConfig[i*3+j].y;
				hAtom[i*3+j].z = refConfig[i*3+j].z;
			}
		}
	}

	//compute the acceptance ratio and new dV
	acceptRatio =(double)acceptCountNPT/(double)moveCountNPT;

	if ((moveCountNPT)%40==0 && simulationType== "EQ") {
		xdl = (acceptRatio - 0.4)/0.4;
		deltaLmax = deltaLmax + xdl*deltaLmax;
	}

	//calculate density
	if (simulationType=="PROD"){
		rho = ((double)nMolecules*44.011)/(Vnew*1000000.0*N_A);
		//printf("%d rho = %lf\n", moveCountNPT, rho);
		sumrho = sumrho +rho;
		avgrho = sumrho/moveCountNPT;
		//	printf("avgwho = %lf\n", avgrho);
	}

	//print out log file
	//	if(moveCountNPT%1 ==0){
	// prec hp = enthalpy/ENERGY_CF;
	nptLog << moveCountNPT << "\t\t" << rho << "\t\t" << avgrho << "\t\t" << potentialEnergy/ENERGY_CF << "\t\t" << acceptRatio << "\t\t" << box.x << "\t\t" << enthalpy << endl; 
//	} 
}

